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The phenomenon of PP (parity- and time-reversal) symmetry breaking is conventionally 
associated with a change in the complex mode spectrum of a non-Hermitian system that 
marks a transition from a purely oscillatory to an exponentially amplified dynamical regime. 
In this work we describe a new type of PP -symmetry breaking, which occurs in the steady- 
state energy distribution of open systems with balanced gain and loss. In particular, we 
show that the combination of nonlinear saturation effects and the presence of thermal or 
quantum noise in actual experiments results in unexpected behavior that differs significantly 
from the usual dynamical picture. We observe additional phases with preserved or ‘weakly’ 
broken PP symmetry, and an unconventional transition from a high-noise thermal state to 
a low-amplitude lasing state with broken symmetry and strongly reduced fluctuations. We 
illustrate these effects here for the specific example of coupled mechanical resonators with 
optically-induced loss and gain, but the described mechanisms will be essential for a general 
understanding of the steady-state properties of actual RT-symmetric systems operated at 
low amplitudes or close to the quantum regime. 
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I. INTRODUCTION 


In 1998 Bender and Boettcher [1] described a class of non-Hermitian ‘Hamiltonians’ that ex¬ 
hibit a purely real energy spectrum, a surprising fact which they attributed to the underlying VT 
(parity- and time-reversal) symmetry. Their observation triggered considerable interest in discrete 
and continuous systems with VT symmetry along with alternative non-Hermitian formulations of 
quantum theory [2]. Such fundamental considerations remain speculative, but there exist many 
classical systems in which "PT-symmetric dynamics can be obtained with appropriately engineered 
loss and gain. This was first pointed out in the context of photonic waveguides [3-6], lattices, and 
resonators [7-10]. Other examples include cold atoms [11-13] and optomechanical devices [14-16]. 
Of particular interest in such systems is the breaking of VT symmetry, i.e., when by tuning a param¬ 
eter the energy spectrum becomes complex and the eigenvectors no longer exhibit the underlying 
VT symmetry. This phenomenon was first experimentally observed in optical waveguides [5, 6], 
and is currently the subject of intense experimental and theoretical research. 


In this work we go beyond this dynamical picture and address an interesting and still open 
question: what are the steady states of actual physical systems with VT symmetry? This ques¬ 
tion becomes especially important for atomic or microscopic solid-state realizations of gain-loss 
systems. Here nonlinear saturation effects as well as the presence of thermal and quantum noise 
have a crucial influence on the system’s dynamics and the long-time behavior can no longer be 
inferred from an eigenvalue analysis only. By focusing on the experimentally relevant example 
of coupled mechanical resonators with optically-induced gain and loss (see Fig. 1) we show that 
■PT-symmetry breaking in the steady state exhibits various unexpected features and in general 
occurs via additional intermediate phases with retained or ‘weakly’ broken VT symmetry. Most 
importantly, we identify an unconventional transition from a high-noise balanced energy distribu¬ 
tion to a parity-broken lasing state with strongly reduced fluctuations. This transition generalizes 
the phenomenon of PT-symmetry breaking—hitherto defined only for eigenstates—to steady-state 
distributions of noisy systems. The mechanisms described here occur in systems of two coupled 
modes as well as in multi-resonator arrays, and will thus be of relevance for a large range of VT 
symmetric systems operated at low amplitudes and close to the quantum regime. 
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FIG. 1. a) Setup of two coupled mechanical resonators with optically-induced gain and loss, b) Scheme for 
engineering mechanical gain or loss via an optically-driven three-level defect. Depending on the detuning 
of the pump (green arrow), phonon-induced transitions between the two near-degenerate excited states jx) 
and \y) lead to a net absorption or emission of phonons of frequency Wm (see reference [ 21 ] for more details). 

II. MODEL 


To motivate the following analysis by a concrete physical system, we consider a setup of two 
micromechanical resonators as shown in figure 1 a). The main findings discussed below, however, 
are more general and can be studied in other equivalent realizations, e.g., with coupled optical or 
microwave resonators. The mechanical resonators have a bare vibrational frequency cUm and they 
are mutually coupled, e.g., mechanically via the support, with strength g. In addition, optical 
or electrical cooling [17-21] and pumping [21-28] schemes are used to induce mechanical loss for 
one resonator and an equal amount of mechanical gain for the other. In a frame rotating with 
Wm, the semiclassical dynamics of the system is then described by the Ito stochastic differential 
equation [29, 30] 

d \ / r+(a) -ig W « ^ / F+{t) \ 

/3/ \ -ig r_(/3 )) \ 13 ) \F.it) ) 

where a and /3 are the dimensionless amplitudes of the pumped and cooled mode respectively 
(more details on the derivation of equation ( 1 ) are given in A). 

The optically-induced gain and loss rates considered here are of the form 


r±(a)= ± 


(1 -b \a\‘^/noy 


■7, 


( 2 ) 


where T is the maximal rate and yTio is the saturation amplitude. The value of v characterizes the 
underlying heating or cooling mechanism and will be treated here as an adjustable parameter. For 
the three-level scheme depicted in hgure 1 b) this parameter takes the value v = 2 [21]. Instead, 
for conventional laser amplification with inverted two-level systems we would obtain v = \ [29]. 
Finally, 7 is the bare mechanical damping rate. Since we are interested in the 'PT-symmetric 
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limit (defined below), we assume 7 /r —)■ 0. However, in all our calculations we retain a finite 
7 > 0 , which describes the actual physical situation and results in a well-defined steady state for 
all parameter regimes. 

In equation (1) the (complex) stochastic forces F±{t) represent two independent white-noise 
processes with {F^{t)F±{t')) = D±6{t — t'). For resonators coupled to a reservoir of temperature T 
the diffusion rates are F>+(a) = Dq{a) -|- 27 A^th and = 27 iVth, where iVth = _ l)“i. 

The contribution Dq(a —)• 0) = 2r for the gain mode represents the intrinsic quantum noise 
associated with any amplification process and suggests that noise is a fundamental property of 
PT-symmetric systems [31-33]. 


III. TT-SYMMETRY BREAKING 


By ignoring the effect of noise and assuming |ap,|/lp <C no, equation (1) can be written as 
dtil^ = —iFlip, with a state vector if; = (a, /3)^ and a non-Hermitian ‘Hamiltonian’ 


H = 




( 3 ) 


This Hamiltonian is invariant under a combined parity {V) and time-reversal (7~) symmetry for all 
values of r/g, where V : {a, -H- {(5,a)'^ simply corresponds to an exchange of the two modes 

and T : i ^ —i. The eigenvalues and (unnormalized) eigenvectors of H are [34] X± = ±\/< 7 ^ — T^, 
■ 0 + = (e*2,e“*2)'^ and V’- = (*e“*2 ^ where sin(0) = T/g. We see that despite being 

non-Hermitian, for T < g both eigenvalues X± are real, indicating purely coherent oscillations 
between the two modes. In this regime of unbroken VT symmetry the eigenvectors are eigenstates 
of the symmetry operator, i.e., 'PF'il’± = 'ip±. For F > gf, both eigenvalues become imaginary and 
correspond to a gain and a loss eigenmode. In this regime 6 is complex and 'il>± no longer possess 
the same symmetry as H. One thus speaks of a "PT-symmetry-breaking transition occurring at 
the (exceptional) point F = g [2] . 

This conventionally studied PT-symmetry-breaking effect captures the change in the system’s 
transient dynamics, as it is observed, e.g., in the propagation of coupled optical modes through 
an active medium with loss and gain [5, 6 ]. However, a simple eigenvalue analysis does not afford 
any conclusions concerning its steady state. Firstly, due to an exponentially amplified mode in 
the symmetry-broken phase, nonlinear effects must be taken into account in order to determine 
the system’s long-time behavior [4, 35-38]. Secondly, due to non-decaying oscillations in the PP- 
symmetric phase, any source of noise would heat up the system indefinitely, and within a linearized 
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FIG. 2. Steady state V'ss = (ass,/3ss)^ of the T^T-symmetric phonon system in the absence of noise. Steady 
state occupation numbers of the two modes for a) = 2 and b) = 1 and 7/3 = 10“^. In the limit-cycle 
phase IIIw both modes oscillate over the range indicated by the shaded area, c) Illustration of the relaxation 
dynamics of \a\^ (red) and |/3p (blue), d) The resulting steady state (green square) for each phase. 

model the steady state would again be ill-defined. 


IV. STEADY-STATE iPT-SYMMETRY BREAKING 

Let us now return to the full model in equation (1) and evaluate the stationary state V’ss = 
{ass, Pss)'^ in the absence of noise. Figures 2 a) and b) show the mode occupation numbers |assp 
and |/3ssP as a function of T/g, and for the two relevant cases = 2 and z^ = 1. Firstly, we 
observe in both plots the expected transition at r/ 5 r|i_>ii = 1, below which (phase I) the system 
dynamics is oscillatory with a small overall damping to |assP = |/3ssP = 0 due to a finite 7 > 0 . 
Above this transition point (phase II) the linearized system dynamics becomes unstable and both 
resonators reach a finite steady-state occupation number |assP/uo = |/3ssP/?^o = — 1 + 

0 ( 7 ), determined by the saturation of r±(a). However, this steady state is still an eigenstate of 
the symmetry operator, VTipss oc tpss, and contrary to our naive expectation the system remains 
"PT-symmetric beyond the conventional transition point. 

The existence of a PT-symmetric steady-state with non-vanishing amplitude is permitted in 
models like the present one, where PP-symmetry is retained even in the nonlinear regime. This 
means that a steady state V’ss with |ass| = |/3ss| p 0 results in an equal suppression of both the 
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gain and the loss rate, r_|_(Q;ss) = —r_(/?ss) [see equation (2)]. This implies that there exists a 
symmetric state which satisfies ^ss = 0 for all values of T /g. However, as one can see in figure 2, 
at larger values of T the system eventually switches to a different, symmetry-broken state. As we 
discuss now, the details of this "PT-symmetry-breaking mechanism and the resulting state depend 
on the actual form of the saturable gain, which is here determined by the parameter v. 

Considering first the case u = 2, we find a second transition sXT/g = A, beyond which (phase 


III) 


as 


l/3s, 


= no X 


' r±y/r{r-Ag) 

2g 


- 1 


( 4 ) 


and the steady state breaks parity, i.e., |ass| / |/3ss|- The PT-symmetry-breaking point can be 
obtained from a linear stability analysis for phase II and for the present model we find 




n- 1 


( 5 ) 


This result would imply the absence of symmetry breaking for n = 1. Instead, for this case we 
observe a Hopf bifurcation around T/g ~ 5.2, beyond which the system approaches a limit cycle 
and undergoes periodic oscillations with the characteristic frequency lOqsc ~ — d) /T. The 

limit cycle is formed symmetrically around a central point V'c with lad = |/3c|, meaning that VT- 
symmetry is preserved on average, but not at each point in time. To distinguish this behavior from 
the previous case we say that in this phase (IIIw) the PT symmetry is only ‘weakly’ broken. In 
contrast to full symmetry breaking, this transition is related to the finite asymmetry induced by 
7 > 0, and would be hidden in the idealized models where 7 = 0. The value for the transition 
point is 

f v + 2u‘^ + \/2j7+~3zP \ 


2 z/2 - 1 


( 6 ) 


y ii->iiiw \ 

independent of the precise value of 7 <C T. The general expressions in equations (5) and ( 6 ) show 
that for all intermediate cases, 1 < < 2 , we obtain 1 < (r/( 7 )ii_>iiiw < (r/( 7 )ii_i.iii < 00 and 

PP-symmetry breaking occurs in three steps I —II —)■ IIIw —)• III. More details on the derivation 
of the above results and the stationary phases for general u are given in B. 

While similar nonlinear phenomena are in general expected for gain-loss systems [16, 37, 39, 40], 
our specihc interest here is to understand the role of dynamical instabilities in the breaking of a 
steady-state symmetry. In particular, the above analysis shows that PP-symmetry breaking occurs 
even in systems where the symmetry is fully retained in the nonlinear regime and a symmetric 
steady-state would be permitted in principle for all parameters. 
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FIG. 3. T^T-symmetry breaking in the limit of large thermal noise, iVth 3> no- a) Steady-state distributions 
of a and P for 7Vth/no = 10 and 7/5 = 10 “^. The values of tq and zg (solid lines) represent the radial 
distance of the distribution maxima from the origin and shaded areas indicate the range of fluctuations, b) 
Plot of the T^T-symmetry parameter A defined in equation ( 8 ). c) Steady-state distribution of a (red dots) 
and (3 (blue dots) in the thermal {T/g = 2, left plot) and in the symmetry-broken {T/g = 10 , right plot) 
phases. 


V. iPT-SYMMETRY BREAKING IN NOISY SYSTEMS 

We now show how the above picture changes in the presence of noise. For clarity we first restrict 
ourselves to a system which is dominated by thermal diffusion, i.e. D± = D = 27Wh- For F = 0 
thermal noise induces additional amplitude fluctuations of about (Aa)^ « 0 /(2^) = and we 
expect that as long as Nth < no the characteristic features shown in hgure 2 , which scale with the 
saturation number no, will only be smeared out, but not change significantly. This is conhrmed 
numerically (not shown) and means that hgure 2 is a good representation of the steady states of the 
system in the weakly nonlinear or low-noise regime. Therefore, we will now address the opposite 
regime Nth > no- 

Figure 3 a) shows the results of a numerical simulation of the stochastic equation (1), from which 
we obtain the steady-state distribution Pss{a,/3) for n = 2 and Nth/no = 10. In the following we 
write a = and /3 = and make use of the fact that the system dynamics is invariant 

under a global phase rotation. The exact marginal distributions F’ss(«) and PssiP) are then htted 
by approximate distributions of the form 

(r — (2 —zpF 

P{a) ^ re ^ P(/3) ~ ze , 


( 7 ) 
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which allows us to extract a radial shift tq and zq and the range of fluctuations (Aa)^ and (A/?)^ 
for both modes (see C). From these values plotted in figure 3 a) we see that the thermal noise now 
completely washes out the features associated with the PT-symmetric phases I and II, and for a 
large range of F the system reaches a steady state (phase T), which is to a good approximation 


thermal, i.e., tq = zq = 0 and (Aa)^ ~ (A/3)^ ~ Nth- Only after a critical value of T/g\ 


T-s>III 


7.5 are the fluctuations suddenly strongly suppressed. In this regime the system relaxes into an 
asymmetric coherent state with rg > zq approximately given by the amplitudes | ass | and | /3ss | given 
in equation (4) and (Aa)^, (A/3)^ ~ yiVthF/f/^, yAth/F <C 1. 

Before we proceed, let us connect this transition to the phenomenon of 'PT-symmetry breaking— 
hitherto defined only for individual states. To do so we introduce the PT-symmetry parameter 


A = 


(l«P)ss + 


< 1 

2 \ — ’ 


( 8 ) 


which vanishes for a random set of states ijji = {ai,/3i)'^ if and only if each state satisfies VTif^i = 
with some real phase Ot- Figure 3 b) shows that indeed A changes at the transition point 
from A ~ Ath = 0.215 for a thermal state to A —)• 1 in the symmetry-broken phase. Note that 
also in the low-noise limit we obtain A = Ath > 0 in phase I and therefore only phase II, where 
A ~ 0, has a strictly PP-symmetric steady state. 

One of the most striking features visible in figures 3 a) and c) is that in sharp contrast to 
a conventional lasing transition, the emerging coherent-state amplitudes after the PT-symmetry 
breaking point are even smaller than the original level of thermal noise. This surprising effect can 
be understood as follows. Although at each instant in time the amplitudes a{t) and /3(t), and 
therefore the gain and loss rates F_|_(a) and F_(/3), can be quite different, the average dissipation 
rate F = (F_ — F_|_)gs when evaluated in the thermal phase vanishes, F ~ 0(7) 0. What 

remains (on average) is the weak coupling to the high-temperature environment. In contrast, 
in the symmetry-broken phase we have (|a|)ss P> (|/3|)ss- Therefore, there is a strong imbalance 
between loss and gain on average, i.e., F ~ 0 {T) > 0, and the resulting net cooling effect suppresses 
fluctuations. Thus, this transition in the average dissipation rate of a stationary system can be 
seen as the counterpart of the transition from real to imaginary eigenvalues in the conventional 
definition of PP-symmetry breaking. What we are still missing, however, is a simple criterion, 
which tells us why the system favors one or the other steady state. 

To clarify the mechanism behind the symmetry-breaking transition we focus on the symmetry- 
broken regime F/^r P> 1, where we can assume that the amplitude of the gain mode |a(t)| ~ |ass| 
and the relative phase ~ 7r/2 are approximately constant. We then obtain the equation of motion 
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FIG. 4. a) Effective potential for the amplitude z = |/3| of the loss mode for T/g = 12. b) Dependence of the 
symmetry breaking point F /g\T^ni on A^th- The solid line represents the prediction from our analytic model 
and the squares show the transition points obtained from the condition A = 0.5 in numerical simulations. 


for the amplitude of the loss mode, z = |/3| [41], 

dtz =-dzU{z) + (9) 

where {r}zit)r]z(t')) = 6{t — t') and (for 7^0) 

U{z) = ~ 2 (i -g|«ssksin((?i>) - ^^^log(z). ( 10 ) 

The function U{z) is an effective potential for z, which is sketched in Fig. 4 a). This potential has a 
local minimum at zq = |/3ss| (corresponding to the steady state given in equation (4) for Nth 0 ), 
which is separated by a finite barrier AU from the unstable region z > Zmax- In the presence 
of noise, a system initially located at z ~ zq can escape over this barrier via thermally activated 

2AU _ 

processes with a characteristic rate i?esc — Rqg , where Rq = —U”{z m\n )U"iz rm^^ )/ ( 47 r^) [30]. 
This rate increases as T is reduced and once i?esc exceeds the bare damping 7 , any configuration 
with fixed a and f3 is rapidly destabilized and the transition to a quasi-thermal state with strongly 
fluctuating amplitudes occurs. In figure 4 b) we compare the transition point r/g'|T-s>iii obtained 
from the condition 7 = i?esc with the numerically evaluated values for various Nth/nQ 1 :^ 1. The 
plot shows that "PT-symmetry breaking in the large-noise regime is very well described by this 
thermal activation model. 


VI. QUANTUM NOISE LIMIT 

As it has been pointed out above, the implementation of any gain mechanism is necessarily 
accompanied by a minimum amount of quantum noise, Pq, which ensures, e.g., the preservation 
of Heisenberg’s uncertainty relations in a quantum mechanical amplification process. For VT- 
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FIG. 5. T^T-symmetry breaking in the quantum-noise limit. The radial position of the steady-state distribu¬ 
tion maxima tq and zq (solid lines) and the range of fluctuations Aa and A/3 are plotted for the case iVth = 0 
and Dq{a) as defined in equation (11). The other parameters for this plot are uq = 10 and j/g = 10“^. 


symmetric systems this implies that noise is not only an intrinsic property, but also that the 
symmetry in the gain and loss modes is broken on a fundamental level. 

To illustrate the effect of pure quantum noise on the system’s steady state we set iVth = 0 and 
consider the quantum diffusion term 


Dq{a) 


2T 

(1 -h laP/no)^’ 


( 11 ) 


only. The assumed form of Dq(a) applies for the three-level gain mechanism shown in figure 1, i.e., 
for the case v = 2 (see A for additional details). Similar to the thermal noise limit discussed above 
we represent the resulting steady-state distributions for the loss and gain modes by the distribution 
maxima rg, zq and the range of fluctuations Aa and A/3 respectively. The results are shown in 
figure 5 for a saturation parameter no = 10. We see that while the steady-state distributions for low 
T do no longer represent thermal states, there still exists a sharp transition between a high-noise 
PT-symmetric phase and a PT-symmetry-broken phase with strongly reduced fluctuations. Thus, 
we conclude that the 'PT-symmetry-breaking mechanism identified above is not very sensitive to 
the details of the noise process and will exist in both thermal and quantum noise limited systems. 
Note that the semi-classical descriptions used in the current analysis is only valid for no ^ 1 while 
for no ~ 1 where quantum effects are most important a full quantum mechanical treatment must 
be employed. Such an analysis is, however, beyond the scope of the current work and will be 
presented elsewhere. 
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FIG. 6. a) Schematic representation of an array of coupled resonators with a T^T-symmetric unit cell, b, 
c) Steady state of a T^T-symmetric array of iV = 12 resonators in b) the low noise regime, A^th/»T-o = 0.01 
and c) the large noise limit iVth/no = 10. In b) the different phases are characterized by the symmetry 
parameter A of a single unit cell and the white-dashed lines indicate the phase boundaries obtained from a 
plane wave ansatz. In c) the solid line shows the analytic prediction for the phase boundary between the 
thermal and symmetry-broken phase and the green squares are the transition points obtained numerically. 
For both plots 7/5 = 10“^ and 1 / = 2, such that the line g' = 0 corresponds to the setting considered in 
figures 2 a) and 3. 

VII. ARRAYS 


Finally, to show that PT-symmetry breaking in the steady state exists also for extended systems, 
we generalize the analysis above to coupled resonator arrays with alternating gain and loss [42, 43], 
as depicted in figure 6 a). In this case the amplitudes and /3n for each unit cell obey 



r+(a„) -ig 
-ig r_(/3n) 



( 12 ) 


where g' is the coupling between the unit cells and Fn,±{t) are independent thermal noise processes. 
Figure 6 summarizes the numerical results for the steady state of an array of V = 12 resonators 
with periodic boundary conditions and u = 2. The observed features can be understood from 
the plane wave ansatz [42, 44] an = /3n = where k = Anj/N. This ansatz maps 

equation (12) onto a two-mode problem for and which is equivalent to equation (1), but with 
the replacement g gk = \g + g'e^^\. This implies that the largest ratio T/\gk\ is always achieved 
for the k = TT mode [43] , which therefore determines the symmetry-breaking properties of the array. 
For the special case g = g', i.e., gk = 0, the gain and the loss modes completely decouple, and 
T^T-symmetry breaking already occurs for F —?■ 0 [9]. For all intermediate parameters the phase 
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boundaries in figure 6 are obtained from the analytic results for the two-mode problem, but with 
g replaced hy g — g'. We see that the single-mode ansatz captures well the relevant physics both 
in the low- and high-noise regime. Note that, however, in the ‘thermal’ phase the behavior of the 
array can actually be much more complicated, since the system may undergo noise-induced jumps 
between multiple metastable configurations. 

VIII. CONCLUSIONS 

We have analyzed the breaking of VT symmetry in the steady state of coupled mechanical 
systems with loss and gain. We have shown how saturation effects and the influence of thermal 
or quantum noise lead to surprising new effects, which are not expected from the conventional 
eigenvalue analysis of idealized "PT-symmetric systems. These findings constitute a significant and 
important first step towards a general understanding of the behavior of realistic systems with VT 
symmetry, especially for systems operated close to the quantum regime or systems with strong 
nonlinearities. 
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Appendix A: Semi-classical laser theory 

In this section we outline the main steps in the derivation of the stochastic differential equa¬ 
tion (1). To do so we first consider an interaction between a single mechanical mode and a second 
auxiliary system, for example an optically pumped defect center [21] or quantum dot [17], which 
is used to engineer mechanical gain or loss. We assume a general coupling of the form 

Hint = A , 


(Al) 
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where a is the bosonic operator of the mechanical mode and A is an operator of the auxiliary 
system. In the absence of the coupling the auxiliary system relaxes with a rate Fa A into 
a steady state described by the density operator p^. This separation of timescales allows us to 
adiabatically eliminate the dynamics of the auxiliary system and derive an effective equation of 
motion for the mechanical mode. Following the standard treatment of semi-classical laser theory 
(see, e.g., section 9.3.2 of reference [45]), the result of such a calculation is a Fokker-Planck equation 
of the form 


P(a, t) = 


d 


da 


d \ 

r(a) + ,^D{a) 


P{a, t). 


(A2) 


da* da 

Here P{a) defined via Pm = f d?a P{a)\a){a\ is the Glauber-Sudarshan phase-space representation 
or P-distribution [29, 45] of the mechanical density operator p^. In this equation the total gain 
function 


F(q;) = Fq(a) - 7, Fq(a) =(A3) 

and the total diffusion function 

L»(a) = 27iVth + Pq(a), Dq{a) = dr (^{A\t)A),, -, (A4) 

are defined in terms of the unperturbed expectation value and correlation function of the operator 
A, where (... )ss = Tra{... Pal- The damping rate 7 and diffusion rate 27Ath encapsulate the effect 
of thermal noise from the environment, whereas Dq(a) is the intrinsic quantum noise associated 
with the gain process. 

Equation (A2) is a standard Fokker-Planck equation and as such it is equivalent to an Ito 
stochastic differential equation with drift F(a) and diffusion D{a) [30, section 4.3.5]. That is, the 
dynamics of the mechanical system may be modeled using the Ito stochastic differential equation 
thus 


da = T{a)adt + \/ D{a)dWc, (A5) 

where dWc = 2~^l‘^{dW\ -|- idW 2 ) is a complex Wiener increment; dWi and dW 2 are independent 
standard Wiener processes [30, section 3.8.5]. By extending this analysis to two modes with mutual 
coupling g we then obtain the two coupled equations 

da = T^{a)adt — igjddt + ZI+(a)dVFc,+, (A6) 

d/3 = F_ {fj)j5dt — igadt + \/ P- {l3)dWc-. (A7) 

By formally introducing the white-noise forces F±{t) = dWc,±/dt, equations (A6) and (A7) repro¬ 
duce equation (1). 










14 


1. Drift and diffusion for a phonon laser induced by a three-level defect 


Reference [21] considers a nitrogen-vacancy (NV) defect centre, which is coupled via strain to 
the bending motion of a diamond nanobeam. The NV center is modeled as a three-level system 
with a ground state \g) and two near-degenerate electronically excited states |x) and \y). If the 
excited-state splitting AExy matches the vibration frequency of the phonon mode Wm then resonant 
transitions are induced between the two states. By optically exciting the upper level \x) this process 
excites the phonon mode (i.e., effects gain) whereas optically exciting the lower level \y) leads to 
a cooling process (see Fig. 1 c) in the main text). In the gain configuration the system can be 
described by the Hamiltonian (h = 1) 

H = UmO^a - AExy\y){y\ + ^ (|fi')(3;| + \x){g\) + A (\y){x\a^ + a|x)(y|^ , (A8) 


where A is the strain coupling constant. For the cooling configuration we obtain the same type of 
coupling, but with the role of \x) and |y) reversed. By assuming a radiative decay rate Fa for both 
of the two excited states, Wm = AExy (gain configuration) and weak driving H <C Fa we obtain [21] 


F 2F 

^ (l + |a|2/no)2’ ^ (l + |a|2/no)3 

where F = /Tl_ and no = Fa/dA^. 


(A9) 


Appendix B: Stability analysis in the absence of noise 

In this section we outline the stability analysis used to obtain the steady-state phases of two 
coupled resonators in the absence of noise, i.e., with D± = 0. By moving to polar coordinates 
a = and /3 = , equation (1) can be rewritten as 

r = —[y — T{r)\r — gs\n{(l))z, (Bl) 

z = -[^+ V{z)\z + gs\n{(i))r, (B2) 

^ = 9 (- - -) cos((/)), (B3) 

where cj) = 6a — 6y- Note that the system is invariant under a combined rotation of a and {3 and 
therefore the evolution of the total phase 6 a + 6p can be neglected. For the last equation we see that 
there are two fixed points for the phase, (/>ss = ±7r/2. Due to finite 7, the stationary occupation 
number of the gain mode \a\‘^ is always slightly larger than that of the loss mode |/3p, wherefore 
'/’ss = 7r/2 is the stable. We therefore set iji = 7r/2 and henceforth study the two-dimensional 
dynamical system with variables r and z. 
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In the following we apply standard dynamical analysis to understand the phases mentioned 
above and predict the transition points. To do so we hrst evaluate the possible stationary solutions 
Tss and Zss of equations (Bl) and (B2), which are given by the solutions of 

gzss = [r(rss) - 7]rss, gvss = [r(2;ss) + (B4) 


The stability of these fixed points is then analyzed using the trace-determinant plane of the Jaco¬ 
bian, i.e. the dynamical matrix of the system linearized about said stationary state. The Jacobian 
for our system is 


J = 


T'{rss)rss + r(rss) - 7 -g 

9 r (2;ss)zss r(zss) 7y 

where the prime denotes the derivative. The trace and determinant of J are 


(B5) 


r = Tr J = -27 r'(rss)rss - r'(2;ss)zss + r(rss) - r(2;ss), (B6) 

J = det J = - [r'(rss)rss + r(rss) - 7][r'(2;ss)^;ss + '['{zss) + 7] (B7) 


respectively. Since the eigenvalues of J may be written entirely in terms of r and 6 thus 

A± = -4J), (B8) 

evaluating r and J at a particular stationary state fully characterizes its stability and local dy¬ 
namical structure. For example, if the real part of both A+ and A_ is negative then the stationary 
state considered is stable, and a non-vanishing imaginary part effects a curl in the phase portrait 
(the r-z-plane). We summarize this characterization in figure 7 a). 


1. Phases 

Phase I: 

The state r^s = ^ss = 0 is an obvious stationary state of equations (Bl) and (B2). Substituting 
this into equations (B6) and (B7) yields (to first order in 7) 

r =-27, 5 = g‘^-V‘^. (B9) 

Consulting figure 7 a), we see that this corresponds to a stable spiral for T < g, but for T > (7 
becomes a saddle node. We thus conclude that phase I is exhibited for 

T 

0 < - < 1 . 

9 


(BIO) 
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FIG. 7. a) Trace-determinant plane of the Jacobian for a two-dimensional dynamical system with cartoons 
of the local phase portraits (reproduced from [46, figure 5.2.8]). The blue- and red-shaded areas indicate 
stable and unstable regions respectively, b) Cartoons of the phase portraits for the four different phases. 
The dashed lines denote nullclines (the curves in the r-z-plane specified by setting either r = 0 or i = 0; the 
intersections of nullclines correspond to stationary states), and the solid lines with arrows denote example 
integral curves (actual trajectories in r and z). c) Cartoon of the general bifurcation diagram for r for 
1 < < 2. Solid blue lines denote stable stationary states, dashed red lines unstable stationary states, and 

the shaded blue and red regions denote stable and unstable limit cycles respectively. This plot has been 
made using data produced by the software MatCont [47] and we have chosen v = 1.5 and 7/5 = 0.001. We 
note the presence of intermediate phases in the transition from phase IIIw to phase III, the green-shaded 
region. In this region the stable stationary state is such that r < z, which vitiates our assumption jaj > j/3j 
that we used to eliminate the phase; for ja] < 1/3] the phase (j) = 7r/2 is unstable. We therefore expect some 
kind of oscillation in this intermediate regime, which is also observed in numerical simulations. 

Phase II: 


In phase II the two occupation numbers are roughly equal, however simply assuming Ugg = Zg, 
yields an inconsistency due to finite 7. Let us therefore substitute the ansatz r^s = + 'yAl'^ +.. 

and Zss = ris'^ J- + ■ ■ ■ and equate alike-orders in 7. The result is 


(Bll) 


-■£’ = X - 1, 

= -xF = +r(r2l) +9|-'. 

We neglect terms of order 7^ and higher. Note that, this stationary state only exists for T/g > 1. 
Substituting equation (Bll) into equations (B6) and (B7) yields 

^ r'(rL°Vis^ + r(r®)+5 

2 i.[(r/ff)V- - l][( 2 z. - l)(r/ff)V- - 2(z. + 1)] 


= -27 1 - 


2 u{T/gy/- + 2{l-u){T/g)V- 


(B12) 

(B13) 


























17 


and 


5 = g^- [r'(rW)r(0) + r(rW)]2 = - g^[ 2 ^ - l - 2v{T/ (B14) 

We see that 5 > 0 only for 1 < F/g' < [ly/iv — 1)]^. However, unlike the previous analysis of 
phase I, in this case r can be positive whilst 6 is positive. Such is the case if — 1)]^ > 

[{u + + 3u'^)j{2u‘^ — 1)]^". Consulting figure 7 a), one sees that as F/g' is increased the 

stationary state (Bll) goes from a stable spiral to an unstable spiral and then to a saddle node. 
We thus conclude that phase II is exhibited for 


F 

1 < — < min 
9 


V + 2v^ + y/2v + 
21^2 _ 1 



(B15) 


Phases IIIw and III; 


In both phase IIIw and III the occupation numbers are not even roughly equal. To analyze this 
regime it is instructive to consider the nullclines (the curves in the r- 2 ;-plane specified by setting 
either r = 0 or i = 0 ; the intersections of nullclines correspond to stationary states): 


z = /(r) - 


7 

-r, 

9 


r = f{z) + -z, 
9 


fir) 


F r 
g (1 + r2/no)'^' 


(B16) 


The function /(r) is always positive, and only zero at r = 0 and r —>■ oo. The maximum is 
located at r/y7io = ( 2 z^ — where /(r) = {T/g){2i')~^{2i' — additional terms 

~il/9)r and {'^lg)z approximately rotate the r- and z-axes for 7 <C g. The possible intersections 
of the two nullclines as a function of T/g are shown in figure 7 c): one intersection (I in the 
figure), two intersections (II and IIIw in the figure), or four intersections (III in the figure). It is 
clear that phase IIIw corresponds to the case of two intersections, and from the previous section 
we know that if [{v + 2v‘^ — yj2u + 3 i/ 2 )/ [2u‘^ — 1)Y < T/g < [v/{v — 1)Y then neither of these 
two intersections corresponds to a stable stationary state. Therefore, by the Poincare-Bendixon 
theorem [46, section 7.3], a limit cycle must exist. We thus conclude that phase IIIw is exhibited 
for 


/ + 2i/2 -)- Y2v P 3i/2 

2i/2 - 1 



(B17) 


but the limit cycle persists slightly beyond this upper bound. Since it is only possible that [i//(z/ — 
1)]^ > [{u + 2 z /2 + Y‘2^v + 3 i/ 2 )/ (2z/^ — 1)]^ if 1 < 1 / < 2, this phase cannot be observed for u >2. 
The limit cycle does not admit a simple analytic form, however by assuming that it is small and 
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FIG. 8. Histograms of the real amplitude r = |a| of the resonator with gain and the fitting distribution (red 
line)—described by equation (Cl) —for three different values of the ratio T/g. a) Fitting distribution in the 
thermal regime, for T/g = 3. b) Fitting distribution in the vicinity of the thermally-activated transition, for 
T/g = 6.3. c) Fitting distribution in the lasing regime, for T/g = 8.6. Note that the y-axis of the last plot 
has been scaled by 1/3. 


centered on stationary state (Bll), which is the only unstable spiral for this regime, one may 
approximate its frequency Wqsc by the imaginary part of the eigenvalue of the Jacobian evaluated 
at equation (Bll). One finds Wosc ~ 2(F — s')/F; this result has been numerically verified for 
the case 1 ^ = 1 (not shown). On the other hand, it is also clear that phase III corresponds to the 
case of four intersections. One may show that there are four intersections only / g > [v/— 
and it appears that (at least for v >2) ol the two extra stationary states one must be an unstable 
spiral and the other a stable spiral. We thus conclude that phase III is exhibited for 


F 

- > 
9 


V 


v-\ 


V 


(B18) 


Note that this phase cannot be observed for v = \. For v = 2 one may easily check that, given 
T/g > — 1)Y = 4, the stable stationary amplitudes are 


= {[F + VF(F - 4<7)]/(2ff) - 1}V2, (B19) 

= {[r - v'r(F - Ag)]/{2g) - 1 }V 2 . (B20) 


Appendix C: Numerical simulations 

For the numerical results presented in Figs. 3 and 4 in the main text we have simulated the 
stochastic equations (A 6 ) and (A7) (and the corresponding extended set of equations for the array) 
using the Euler-Maruyama method (see section 15 of [30]). Since we are interested in the steady 
state of the system, we collect data (complex numbers a{t) and /3(t)) after a time of 5 x tq where 
To = 1/7 is approximately the time scale in which the steady state is approached. Over a period of 
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45 X To we select 4000 random data points. We repeat the procedure 80 times with random initial 
conditions. 

From the numerical data we obtain the marginal steady-state distributions Pssict) and Pssif^) for 
the two modes. Since the system is invariant under a combined rotation of a and /3 in phase space, 
also the marginal distributions are radially symmetric. Figure 8 shows examples of the resulting 
radial distribution Pss(r = |a|) for the gain mode, before, close to, and after the transition. 

To obtain a simple characterization of the system’s steady state, we fit the radial distribution 
Pss('i~ = |ct|) by a distribution of the form 

Pro,a{r) = J\f X r X e (Cl) 

where M is a normalization constant. It corresponds to a thermal state for tq = 0 and a = -^/Wh 
and approaches the distribution of a coherent state with random phase in the limit a <C tq. The 
optimal values for tq and a are obtained by minimization of the squared difference 

N 

x\ro,a) = Y,iH^-ProAr^))\ (C2) 

i=l 

where Hi represents the hight of the i’th bar of the histogram. 

While this fitting procedure yields accurate values for the radial displacement of the distribution 
maximum, ro, we find that the corresponding values for the width a do not very well capture the 
broad thermal background in the vicinity of the transition [see figure 8 b)]. Therefore, we use 
instead the quantity 

{^af = {\aA-rl (C3) 

where the average (jap) is calculated directly from the data set, to represent the range of fluctua¬ 
tions. 
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